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Resumen 


Se presenta una metodología para calcular de forma numérica el flujo con 
potencial, empleando mallas curvilíneas compuestas. Las mallas curvilíneas 
tienen la ventaja de adaptarse a las formas irregulares del dominio físico; sin 
embargo, ello no es suficiente en geometrías complejas, por lo cual es necesario 
utilizar mallas compuestas; tal tipo de mallas están conformadas por múltiples 
bloques interconectados. En el presente trabajo se muestra que las técnicas de 
descomposición de dominio proporcionan una metodología eficaz para resolver 
este tipo problemas; se obtiene primero la solución en la interfaz de las mallas 
interconectadas, o bloques, y después la solución en cada una de las mallas. Con el 
fin de validar la metodología, se realiza una comparación con soluciones analíticas. 
Se presenta un caso de aplicación. 


Palabras clave: flujo con potencial, mallas compuestas, mallas curvilíneas, 
descomposición de dominio. 


Introducción cas: numéricas, analíticas, mapeo conforme o 
por superposición de flujos simples, entre otras 
La teoría del flujo con potencial permite (Grenger, 1995). Las técnicas no numéricas 
analizar flujos reales bajo la hipótesis de que 


son irrotacionales. Un flujo real está sujeto a 


tienen la desventaja de aplicarse sólo en algu- 
nos tipos de geometrías. En este trabajo se pro- 
deformaciones angulares y en principio no pone una metodología de solución numérica 
puede ser modelado con la teoría del flujo con 


potencial. Sin embargo, en ciertos casos puede 


que utiliza mallas curvilíneas compuestas, las 
cuales son de uso más general, ya que permiten 
considerarse como irrotacional; por ejemplo, si dar solución en situaciones donde la geometría 
las fuerzas que generan el movimiento (como del dominio es compleja, que difícilmente se 
la gravedad o un gradiente de presiones) son 


grandes, en comparación con las de origen 


resuelven con técnicas no numéricas y que es 
complicado discretizar con mallas de un único 


viscoso (Sotelo, 2008); en este caso, el flujo que 
se encuentra fuera de la capa viscosa se puede 
considerar irrotacional; a lo anterior se suma 
que no existan zonas de separación. Por tanto, 
el análisis de flujo con potencial proporciona 
una buena aproximación a la distribución de 
velocidades del flujo real en muchos casos 
prácticos. 

La solución se obtiene con la ecuación de 
Laplace, que se resuelve con diferentes técni- 


bloque (Thompson et al., 1999). Las mallas 
curvilíneas compuestas consisten en múltiples 
mallas interconectadas; conservan la ventaja 
de que cada bloque es estructurado. Después 
de un estudio bibliográfico se encontró que las 
técnicas de descomposición de dominio (ver 
por ejemplo Toselli y Widlund, 2005) plantean 
una metodología eficaz para dar solución a 
la ecuación de Laplace en geometrías que se 
dividen en dos o más subdominios. 
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Flujo con potencial 


El flujo con potencial permite calcular el 
campo de velocidades a partir del potencial 
de velocidad, denotado como q. El campo de 
velocidades se obtiene con la relación: 


U=-V6 (1) 


En (1), U es el campo de velocidades, 
definido como U=ui +0 . Al combinar (1) y 
la ecuación de continuidad (Granger, 1995) se 
obtiene: 


2 2 
9%, 00 (2) 
dx dy? 


La ecuación (2) es la ecuación de Laplace 
para el potencial de velocidad. Otra manera 
de analizar el flujo con potencial es la función 
de corriente, y, la cual se relaciona con q a 
través de las condiciones de Cauchy-Riemann 


(Grenger, 1995). Este parámetro satisface 
también la ecuación de Laplace: 
2. 2 
E (3) 
dx” dy 


Condiciones de frontera 


Las ecuaciones (2) y (3) están sujetas a condi- 
ciones de frontera tipo Newmann, Dirichlet 
o Cauchy (Saad, 2003). En este trabajo se pre- 
senta el caso que es una combinación de 
tipo Dirichlet en una región del contorno, 
y Newmann en la parte restante, debido a 
que corresponden con los tipos de fronteras 
manejadas (paredes impermeables y fronteras 
abiertas). Cuando se hace un análisis con d, lo 
común es conocer la magnitud del potencial 
en la entrada y la salida, y se manejan como 
fronteras Dirichlet. La parte restante de los 
contornos son paredes impermeables y se 
cumple la condición 04/41 = 0, que es tipo 
Newmamn. Cuando se resuelve con y, se 
conoce el valor de la función de corriente en las 


paredes, por lo que son tipo Dirichlet, mientras 
las fronteras abiertas son Newmanmn. 


Solución empleando coordenadas 
curvilíneas 


En las geometrías con forma irregular es común 
utilizar dos metodologías de solución numé- 
rica: una es con mallas no estructuradas y 
elemento finito (Segerlind, 1984); otra es con 
mallas curvilíneas (Thompson et al., 1999), con 
las que se pueden utilizar diferencias finitas. 
En este trabajo se utiliza la segunda metodolo- 
gía. El empleo de un sistema coordenado 
curvilíneo consiste en establecer una relación 
única entre el espacio cartesiano x — y, y el 
curvilíneo £ — n, de tal forma que un dominio 
irregular en x — y se represente como un área 
rectangular en £ — n. Lo anterior se indica en la 
figura 1, donde el área delimitada por cuatro 
fronteras irregulares en el lado izquierdo tiene 
una representación rectangular en el derecho. 
En la figura 1 se observa que el uso de este 
tipo de mallas (en su forma tradicional de un 
bloque) tiene la restricción de que el dominio 
en el plano cartesiano debe estar delimitado 
por cuatro fronteras, las cuales corresponden a 
los cuatro lados del rectángulo representado 
en el plano curvilíneo; esta restricción se 
elimina cuando se emplean mallas compuestas. 

El sistema curvilíneo requiere de la trans- 
formación de las ecuaciones del flujo para 
referirlas a este sistema. Con este enfoque 
resulta sencillo implementar un esquema de 
diferencias finitas completamente idéntico para 
todos los nodos interiores de la malla debido a 
que ésta es rectangular. La transformación de 
las ecuaciones diferenciales parte de la regla 
de la cadena. Si se considera una función F 
cualquiera, la derivada respecto a x se calcula 
como: 


dx 0 dx on ox 


La derivada respecto a y se calcula de mane- 
ra semejante. Para evaluar estas derivadas se 
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Figura 1. Relación entre el plano cartesiano 
y el curvilíneo. 


requiere conocer d(€,n)/0(x,y), que se obtienen 
a partir del sistema coordenado curvilíneo. La 
base covariante del sistema curvilíneo se define 
como: 


5% 0x2» 0Y 2 
A (5a) 
se 0x2 Oy 2 
a, mn E (5b) 


Los vectores anteriores son tangenciales a 
los ejes curvilíneos; la base contravariante, que 
es normal a los ejes, está dada por: 


E Es De 

1-VY o O 6 
6 y (6a) 

A LES on (6b) 


La base covariante se calcula numérica- 
mente una vez que la malla ha sido generada, 
mientras que la base contravariante se obtiene 
a partir de la primera, con las expresiones 
(Thompson et al., 1985): 


=)(á,xk) (Za) 


2=)7(a,xk) (7b) 


En (7), J es el Jacobiano de la transformación 
y k es el vector unitario normal al plano 
coordenado. Los elementos métricos, necesa- 
rios para la transformación de la ecuación de 
Laplace, se obtienen a través del producto es- 
calar de los vectores de la base contravariante: 


a = '.¿) con i,¡=12 (8) 


En Thompson et al. (1985) se explica la 
transformación de la ecuación de Laplace de 
coordenadas cartesianas a curvilíneas, que 
parte de las relaciones como la indicada en (4). 
Las ecuaciones (2) y (3) quedan expresadas 
en su forma curvilínea no conservativa, 
respectivamente como: 


Ve =4" 9% +2942 9% +42 9% 


dE* don an” (9) 
e) rm) AL o 
0 on 
y 
Vay=4" y ¡E y qe y 2 NN 
2 a 0 2 
0 ¿on on (10) 


e) 2H) o 


Las ecuaciones (9) y (10) son resueltas con 
un esquema de diferencias finitas centradas 
(Leveque, 2007). Teniendo en cuenta que en el 
plano curvilíneo A =1 y An = 1, y sustituyendo 
el esquema de diferencias, (9) se expresa como 
sigue: 


C.0, +C 40, +C,0, +C0, +C.b. + 


+Cl0re HO 0. Hb) =0 sl 
Donde: 

C,=a"+V%/2 (12a) 

C,= a" -VE/2 (12b) 
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C,=4%+V*9/2 (12c) 
C,=a*-Vn/2 (12d) 
C.=-2(a" +47) (12e) 

C=a*l2 (12f) 


Los subíndices de 4 en (11) se refieren a la 
ubicación de los nodos en la malla curvilínea, 
tal como se indica en la figura 2. Al aplicar 
(11) en cada uno de los nodos de la malla 
curvilínea, incluyendo las condiciones de 
frontera, se obtiene un sistema de la forma A4 
= b, cuya solución es el potencial de velocidad. 
Se procede de manera análoga para la función 
de corriente y. 

Con el potencial de velocidad conocido, el 
campo de velocidades se calcula con la ecuación 
(1), transformada a coordenadas curvilíneas: 


Metodología 
Las mallas compuestas consisten en una 


discretización del dominio físico, empleando 
dos o más mallas interconectadas entre sí. 


01 


Figura 2. Plantilla de cálculo. 


Como ejemplo, en la figura 3 se muestra 
un dominio discretizado con una malla 
conformada por nueve bloques. 

El uso de las mallas compuestas requiere 
una metodología para hacer compatible 
la solución en la interfaz de los bloques. 
En este trabajo se emplean las técnicas de 
descomposición de dominio, como se explica 
en Chan y Mathew (1994), y Toselli y Widlund 
(2005), entre otras referencias. 


Técnicas de descomposición de dominio 


Estas técnicas se dividen en dos categorías: 
las primeras son denominadas técnicas con 
traslape, y las segundas, sin traslape. Las téc- 
nicas con traslape se aplican para resolver los 
sistemas de ecuaciones generados al segmentar 
un dominio (2 en múltiples subdominios Q, 
traslapados entre sí; mientras que las técnicas 
sin traslape se emplean cuando el dominio se 
divide en múltiples subdominios, cuya única 
conexión es a través de una interfaz I', tal como 
se indica en la figura 4. 

En este trabajo se propone emplear la 
metodología sin traslape, ya que facilita la 
generación de la malla curvilínea y el manejo 
de la información. Si se considera que, por 
ejemplo, el dominio se descompone en dos 
subdominios, como se indica en la figura 4, al 
utilizar un esquema de diferencias finitas se 
obtiene un conjunto de nodos que pertenecen 


E 


Figura 3. Descomposición de la bifurcación de un río 


empleando nueve bloques. 
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Figura 4. Descomposición de un dominio sin traslape. 


al subdominio (,, otro conjunto de nodos que 
pertenecen a Q, y finalmente un conjunto de 
nodos que pertenecen a la interfaz TI. El siste- 
ma de ecuaciones obtenido tiene la estructura 
siguiente: 


A Ar % E 1 
2 Azr 0) SS Í 2 (14) 
Ar Arz Ar |[0r fr 


Los componentes involucrados en (14) 
son vectores y matrices; A, y A, son las 
matrices resultantes de las discretizaciones 
en los subdominios (2, y Q,, mientras que la 
matriz A, corresponde a la discretización de la 
interfaz IP”. Las matrices A, , y A,, establecen la 
interacción entre la interfaz y los subdominios. 
Finalmente, los vectores q, 4, y 6, son la 
solución en Q,, Q, y I, respectivamente. 

El sistema (14) está conformado por una 
matriz en bloque o matriz de matrices; puede 
resolverse con la eliminación de Gauss por 
bloques (Ortega, 1988), por lo que se obtiene el 
sistema equivalente: 


A; Ar 0 f 
A, Ar 0, > E (15) 
s |lór| 13 


En el sistema (15), la matriz S y el vector g 
están definidos de la siguiente manera: 


S= Ap Ar, Aj Ar Aro AA (16) 


2 =AA Buds ta (17) 
La solución en los nodos de la interfaz es 
obtenida al resolver el sistema de ecuaciones 
del último renglón en el sistema (15): 
SOp=8 (18) 
Una vez que se conoce 0, se obtienen las 
soluciones $, y 6, al resolver los sistemas de 
ecuaciones asociados con los subdominios 
Q, y O, ubicados en los renglones 1 y 2, 
respectivamente, del sistema (15): 


Aj, = (A 5 Ajróp) (19) 


Asd, = (A o Azrbp) (20) 

Con los sistemas (18) al (20) se obtiene la 
solución del dominio completo. Las ecuaciones 
(16) y (17) corresponden a una descomposición 
en dos subdominios; sin embargo, la 
generalización para n subdominios es de la 


siguiente forma: 


(21) 
(22) 


En (21) y (2), que proporcionan la 
matriz y el vector necesarios para obtener la 
solución en la interfaz por medio de (18), se 
observa que es necesario invertir las matrices 
de cada uno de los subdominios. Aunque 
en la literatura (Chan y Mathew, 1994) se 
recomienda resolver (14) de forma iterativa, 
en este trabajo se propone utilizar un método 
de solución directa. Esto involucra invertir las 
matrices de los subdominios para obtener la 
matriz S y el vector g. Esto se debe a que en 
pruebas realizadas en diferentes mallas con 
un número de nodos de algunos miles, no 
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resultó computacionalmente caro trabajar de 
esta manera. Las matrices inversas A,* pueden 
calcularse a través de la factorización LU, 
obteniéndose la matriz inversa de la forma: 


ART=UuTL” 


1 


(23) 


Invertir las matrices A, es siempre posible 
cuando sean no singulares (Saad, 2003), lo 
cual ocurre prácticamente en todos los casos 
para la ecuación de Laplace: los coeficientes 
(12) producen un sistema en banda donde la 
diagonal principal es dominante. Una ventaja 
de utilizar la factorización LU es que, una vez 
que se obtiene la solución en la interfaz a través 
de resolver (18), puede utilizarse el resultado 
de la factorización para resolver el sistema 
siguiente: 

Af; =(f = A,rór) (24) 

El sistema (24) es una generalización de (19) 
y (20). 

Para ejemplificar el proceso de discretiza- 
ción utilizando una malla compuesta, se 
analiza el caso indicado en la figura 5, 
correspondiente a un dominio que se ha 
descompuesto en tres bloques. Los nodos se 
numeran consecutivamente, incluyendo los de 
la interfaz. Con esta numeración se obtiene el 
sistema de ecuaciones: 


A, Aur % f 
A, 0 Ar d, 59 f, (25) 
Az A>r 0d, f, 


Ari Aro Ars r JLór hi Tr 


Además de la numeración global, cada 
bloque y la interfaz tienen numeración local, 
como se indica en el cuadro 1. 

Las matrices componentes de (25) se ob- 
tienen a partir de la información mostrada en 
la figura 5, donde los bloques se representan 
en el espacio curvilíneo; en este ejemplo, por 
simplicidad se tiene como base contravariante 
W=i y d=j, por lo que V?%E=V*N=0, y, 


Figura 5. Dominio descompuesto en tres bloques. 


Cuadro 1. Relación de nodos globales y locales. 


Bloque Nodo global Nodo local 


1 


pues 


Nm 
NV|M|I|D|O|RA2[|o0|n 


Interfase 14 


ejo|n|R|2ejoOo|N|R|2ajo|[N|R|Baj|oj|mn 


como es costumbre en la malla curvilínea, se 
tiene A£ =An=1. Al aplicar los coeficientes 
definidos por (12), para el nodo 1 se tiene que 
0, =0,, 0, =0, y 0, = 0,, que sustituidos en (11) 
da por resultado: 


0, +0, 40,+0,,+0,=0 (26) 
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Los valores cruzados en (11), como por 
ejemplo dé, , se cancelaron en (26) debido a que 
la malla es ortogonal. Las matrices para los 
bloques (25) quedan: 


4 1 1-0 
1. -4 0 
e E 0 -4 1 (27) 
y 11 
=2-.1 0 
ap E AE OO da) 
Ela. 4 <t.1 
0 0 1 -4 
000.0 
O y A 
Ajr= 29a 
GS Go MS 1228) 
01.01 
0010 
a ic (29b) 
“(0.0.0 1 
000.0 
1000 
" 0100 dé 
= Cc 
*- 10-000 
000.0 


En el cuadro 1 se observa que cada bloque 
y la interfaz están formados por cuatro nodos, 
por ello las matrices A, y A, son en este caso de 
4 x 4. Una vez que se tienen las matrices como 
se indican en (27) a (29), se procede a calcular 
las matrices S y g (este último vector se obtiene 
con los componentes f, que se determinan 
con las condiciones de frontera) para obtener 
la solución en la interfaz, y posteriormente en 
cada uno de los bloques. 


Validación del esquema 


Se realiza la comparación de esta metodología 
con resultados analíticos. El caso corresponde a 
la superposición de un flujo uniforme más una 
fuente; la función de corriente es y = Y... 

uniforme 
+ Wien (Granger, 1995), que reproduce el flujo 
alrededor de un obstáculo, como el de la figura 
6a. La ecuación de la función de corriente, en 


coordenadas polares (r,0), es la siguiente: 


y= Ursen(0)+ Lg (30) 
21 


donde U es la velocidad del flujo uniforme y q 
es el caudal unitario en la fuente. Con (30) se 
obtienen las líneas de corriente indicadas en la 
figura 6a, donde se utiliza U =2 m/s y q = 20 
m?/s. El dominio empleado para la solución 
numérica es el delimitado por las líneas de 
corriente y =-190 m?/ s en el sur, y = 210 m?/s 
en el norte, y por las rectas x = 140 m en el este, 
y x = -140 m en el oeste. En el norte y sur se 
tienen fronteras tipo Dirichlet, mientras que 
en el este y oeste son tipo Newman. La parte 
correspondiente al obstáculo se maneja como 
Dirichlet con y = 10 m?/s. Se emplean cuatro 
bloques para discretizar el dominio, tal como 
se indica en la figura 6b. 

Se compara la solución analítica con la 
obtenida numéricamente en la recta A-A', 
indicada en la figura 6a. Los resultados 
analítico y numérico se muestran en la figura 7; 
en términos prácticos, son los mismos. 


Aplicación 


Se estudia la bifurcación del río Mezcalapa 
en los ríos Samaria y Carrizal; en este caso 
se calculó el campo de velocidades y la 
distribución de gastos en los dos ríos de 
salida. La malla se indica en figura 8 y está 
conformada por 1 040 nodos. En la entrada 
del río Mezcalapa el potencial de velocidad 
se estableció como q = 4 000 m?/s, y en las 
salidas de los ríos Samaria y Carrizal como q = 
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—_ _______________________JJJ_ ____ y = -110— 
=-150= 
=-190— 
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Figura 6. a) Líneas de corriente, b) malla compuesta 
de cuatro bloques. 


Función de corriente 


Analítica — 


Numérica 


Línea A-A” 


Figura 7. Comparación de los resultados analítico 


y numérico. 


0 m?/s. Por otra parte, las márgenes de los ríos 
se manejaron como fronteras tipo Newmann, 
esto es, d0/ dn =0. 

Se resolvió para el potencial de velocidad 
debido a que puede asociarse intuitivamente 
con el tirante, h, ya que el sentido del campo 
de velocidades ocurre hacia donde hay un 
gradiente negativo de h. 

La distribución de gastos se muestra en el 
cuadro 2. El campo de velocidades calculado 
se indica en la figura 9. Debe tenerse en cuenta 
que el caudal obtenido es unitario respecto a la 
dirección perpendicular al plano. 


y 
$ 
S 
ES 
A 
9 


S E 


x= 


Figura 8. Malla conformada por nueve bloques de la 
bifurcación del río Mezcalapa. 


Cuadro 2. Distribución del caudal en los ríos. 


Río O (n7/s) % 
Mezcalapa 914 100 
Carrizal 294 32 
Samaria 620 68 


Figura 9. Campo de velocidades. 
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Conclusiones 


La teoría del flujo con potencial sigue siendo 
una herramienta útil cuando la física del 
problema permite considerar al flujo como 
irrotacional. En este trabajo se presenta una 
metodología de solución alternativa cuando 
se tienen geometrías donde es complicado 
utilizar una solución no numérica o que resulta 
difícil de discretizar con una malla de un único 
bloque. 

Un aspecto importante a considerar es que 
aunque la metodología aquí presentada se 
aplicó al flujo con potencial, puede utilizarse 
en cualquier otro tipo de problemas, donde 
el dominio físico requiera una discretización 
basada en mallas compuestas, debido a que 
las técnicas de descomposición de dominio 
permiten dividir un problema grande en 
múltiples problemas más pequeños. Cualquier 
esquema numérico de solución de ecuaciones 
diferenciales parciales que dé como resultado 
un sistema de ecuaciones lineales puede 
resolverse eficientemente con el planteamiento 
de mallas compuestas presentado en este 
trabajo. 
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Mendoza-Reséndiz y Berezowsky-Verduzco, Solución numérica de flujo con potencial en geometrías complejas 


Abstract 


MENDOZA-RESÉNDIZ, A. € BEREZOWSKY-VERDUZCO. M. Numerical solution of 
potential flow in complex geometries. Water Technology and Sciences, formerly Hydraulic 
engineering in Mexico (in Spanish). Vol. II, No. 2, April-June, 2011, pp. 183-192. 


A methodology for the numerical computation of potential flow using composite curvilinear 
grids is presented. Curvilinear grids can be adapted to physical domains with irregular 
geometry. However, in complex geometries this is not enough; in such cases, composite grids 
should be used, which are formed by multiple interconnected blocks. This paper shows that 
the use of domain decomposition techniques provides an efficient methodology for solving this 
kind of problems, obtaining the solution first at the interface, and later in each block. In order 
to validate the proposed methodology a comparison with analytical solutions is performed. An 
application case is also presented. 


Keywords: potential flow, composite grids, curvilinear grids, domain decomposition. 
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